# Load data

d = read.csv((here("final_analyses", "output", "data", # file path
                   "respiratory-data.csv")),           # file name
             header = T)
d = d %>% 
  as_tibble() %>% 
  select(!starts_with("X")) %>% 
  
  # To create columns with the difference between sample size and events per arm
  # This will be useful while using Beta distributions,
  # diff = the Beta (b) parameter in a Beta distribution (a,b)
  mutate(diff_toci = total_toci - events_toci,
         diff_control = total_control - events_control)

Let’s check again how many trials (other than RECOVERY) we have for each outcome and subgroup

d %>%
  filter(!str_detect(oxygen, "pooled"),
         trial != "RECOVERY") %>% 
  group_by(outcome, oxygen) %>%
  arrange(oxygen) %>%
  count() %>%
  rename("number of trials" = n) %>% 
  flextable() %>%
  autofit()
### To create separate objects per subgroup
# Here I rename columns so these variables are usable in my beta_fun()

# Isolating RECOVERY
recovery_simple_oxygen = 
  d %>% 
  filter(trial == "RECOVERY",
         oxygen == "simple oxygen") %>% 
  
  # Only RECOVERY, so there is no need to sum, just to rename the columns
  rename(sum_events_toci = "events_toci",
         sum_total_toci = "total_toci",
         sum_diff_toci = "diff_toci",
         sum_events_control = "events_control",
         sum_total_control = "total_control",
         sum_diff_control = "diff_control")

recovery_non_invasive = 
  d %>% 
  filter(trial == "RECOVERY",
         oxygen == "non-invasive ventilation") %>% 
  
  # Only RECOVERY, so there is no need to sum, just to rename the columns
  rename(sum_events_toci = "events_toci",
         sum_total_toci = "total_toci",
         sum_diff_toci = "diff_toci",
         sum_events_control = "events_control",
         sum_total_control = "total_control",
         sum_diff_control = "diff_control") 

recovery_invasive = 
  d %>% 
  filter(trial == "RECOVERY",
         oxygen == "invasive ventilation") %>% 
  
  # Only RECOVERY, so there is no need to sum, just to rename the columns
  rename(sum_events_toci = "events_toci",
         sum_total_toci = "total_toci",
         sum_diff_toci = "diff_toci",
         sum_events_control = "events_control",
         sum_total_control = "total_control",
         sum_diff_control = "diff_control")

## Group all other trials per outcome and sum the events/sample size

trials_simple_oxygen =
  d %>% 
  filter(trial != "RECOVERY",
         oxygen == "simple oxygen") %>% 
  
  # Only COVACTA, so no need to sum, just to rename the columns
  rename(sum_events_toci = "events_toci",
         sum_total_toci = "total_toci",
         sum_diff_toci = "diff_toci",
         sum_events_control = "events_control",
         sum_total_control = "total_control",
         sum_diff_control = "diff_control")

trials_non_invasive = 
  d %>% 
  filter(trial != "RECOVERY",
         oxygen == "non-invasive ventilation") %>%
  group_by(outcome) %>% 
  summarise(sum_events_toci = sum(events_toci),
            sum_total_toci = sum(total_toci),
            sum_diff_toci = sum(diff_toci),
            sum_events_control = sum(events_control),
            sum_total_control = sum(total_control),
            sum_diff_control = sum(diff_control))

trials_invasive = 
  d %>% 
  filter(trial != "RECOVERY",
         oxygen == "invasive ventilation") %>%
  group_by(outcome) %>% 
  summarise(sum_events_toci = sum(events_toci),
            sum_total_toci = sum(total_toci),
            sum_diff_toci = sum(diff_toci),
            sum_events_control = sum(events_control),
            sum_total_control = sum(total_control),
            sum_diff_control = sum(diff_control))

Mortality

Regarding this outcome, negative risk differences mean benefit.


Using a non-informative prior

Simple oxygen only

### Generate posterior beta distributions using a noninformative beta(1,1) prior

set.seed = 123 # set seed for reproducibility (rbeta() function)

n = 10e4 # To determine sampling size

# Here I use the uninformative_posterior_beta() in which you can find in the
# following path: 
# here("final_analyses", "script", "functions", "beta_distributions.R")

beta_recovery_simple_oxygen = uniformative_posterior_beta( 
  recovery_simple_oxygen, # data
  "mortality" # outcome
) 

### Plot!

# I will use my own functions to standardize my plots

risk_comparison_plot(
  beta_recovery_simple_oxygen, # Data object

  ### start using quotes
  
  "gray50", # Color code for control group
  "#D3A464", # Color code for other group
  "\nRisk (%)", # X axis label
  "RECOVERY trial", # Title
  "Simple oxygen only subgroup\n", # Subtitle
  
  ### stop using quotes

  10, # X axis inferior limit
  30, # X axis superior limit
  2 # X axis spacing for ticks
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlabel = "\nRisk difference (%)"
xlim_inf = -10
xlim_sup = 5

p1 = risk_difference_plot(
  beta_recovery_simple_oxygen, # Data object

  ## start using quotes
  "#7EBAC2", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # X axis label
) 

p2 = cumulative_risk_difference_plot(
  beta_recovery_simple_oxygen, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  1, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "RECOVERY trial",
  subtitle = "Risk difference between groups on simple oxygen only"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

Non-invasive ventilation

### Generate posterior beta distributions using a noninformative beta(1,1) prior

set.seed = 123 # set seed for reproducibility (rbeta() function)

n = 10e4 # To determine sampling size

# Here I use the uninformative_posterior_beta() in which you can find in the
# following path: 
# here("final_analyses", "script", "functions", "beta_distributions.R")

beta_recovery_non_invasive = uniformative_posterior_beta( 
  recovery_non_invasive, # data
  "mortality" # outcome
) 

### Plot!

# I will be using my own functions to standardize my plots

risk_comparison_plot(
  beta_recovery_non_invasive, # Data object

  ### start using quotes
  
  "gray50", # Color code for control group
  "#D3A464", # Color code for other group
  "\nRisk (%)", # X axis label
  "RECOVERY trial", # Title
  "Non-invasive ventilation subgroup\n", # Subtitle
  
  ### stop using quotes

  28, # X axis inferior limit
  48, # X axis superior limit
  2 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlabel = "\nRisk difference (%)"
xlim_inf = -12
xlim_sup = 5

p1 = risk_difference_plot(
  beta_recovery_non_invasive, # Data object

  ## start using quotes
  "#7EBAC2", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # Spacing between ticks in X axis
) 

p2 = cumulative_risk_difference_plot(
  beta_recovery_non_invasive, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  1, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "RECOVERY trial",
  subtitle = "Risk difference between groups on non-invasive ventilation"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

Invasive mechanical ventilation

### Generate posterior beta distributions using a noninformative beta(1,1) prior

set.seed = 123 # set seed for reproducibility (rbeta() function)

n = 10e4 # To determine sampling size

# Here I use the uninformative_posterior_beta() in which you can find in the
# following path: 
# here("final_analyses", "script", "functions", "beta_distributions.R")

beta_recovery_invasive = uniformative_posterior_beta( 
  recovery_invasive, # data
  "mortality" # outcome
) 

### Plot!

# I will be using my own functions to standardize my plots

risk_comparison_plot(
  beta_recovery_invasive, # Data object

  ### start using quotes
  
  "gray50", # Color code for control group
  "#D3A464", # Color code for other group
  "\nRisk (%)", # X axis label
  "RECOVERY trial", # Title
  "Invasive mechanical ventilation subgroup\n", # Subtitle
  
  ### stop using quotes

  36, # X axis inferior limit
  58, # X axis superior limit
  2 # Spacing between ticks in X axis
) + 
  
  # Extra function to better limits the axis in this case
  coord_cartesian(c(36, 58))
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlabel = "\nRisk difference (%)"
xlim_inf = -16
xlim_sup = 12

p1 = risk_difference_plot(
  beta_recovery_invasive, # Data object

  ## start using quotes
  "#7EBAC2", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  4 # Spacing between ticks in X axis
) 

p2 = cumulative_risk_difference_plot(
  beta_recovery_invasive, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "RECOVERY trial",
  subtitle = "Risk difference between groups on invasive mechanical ventilation"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

Using evidence-based priors

Simple oxygen only

priors_simple_oxygen_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "simple oxygen"
  )

priors_simple_oxygen =
  rma(
    # Outcome
    measure = "RD", # risk difference
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = priors_simple_oxygen_raw %>% filter(outcome == "mortality"),

    slab = paste(trial),
    method = "REML"
  )

forest(priors_simple_oxygen)

set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size

xlabel = "\nRisk difference (%)"

data_prior_posterior_plot(
  recovery_simple_oxygen_normal, # Output from normal_approximation()
  
  # start using quotes
           "#364D55", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY and previous evidence", # Title
           "Simple oxygen only subgroup", # Subtitle
           # stop using quotes
           -10, # X axis inferior limit
           20, # X axis superior limit
           5 # Spacing between ticks in X axis
  )
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlim_inf = -8
xlim_sup = 4

p1 = posterior_difference_plot(
  recovery_simple_oxygen_normal, # Data object

  ## start using quotes
  "#A65041", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # Spacing between ticks in X axis
) 

p2 = posterior_cumulative_plot(
  recovery_simple_oxygen_normal, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  1, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial and previous evidence",
  subtitle = "Simple oxygen only subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

Non-invasive ventilation

priors_noninvasive_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "non-invasive ventilation"
  )

priors_noninvasive =
  rma(
    # Outcome
    measure = "RD", # risk difference
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = priors_noninvasive_raw %>% filter(outcome == "mortality"),

    slab = paste(trial),
    method = "REML"
  )

forest(priors_noninvasive)

set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size

data_prior_posterior_plot(
  recovery_noninvasive_normal, # Output from normal_approximation()
  
  # start using quotes
           "#364D55", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution - RECOVERY and previous evidence", # Title
           "Non-invasive ventilation subgroup", # Subtitle
           # stop using quotes
           -15, # X axis inferior limit
           5, # X axis superior limit
           5 # Spacing between ticks in X axis
  )
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlim_inf = -12
xlim_sup = 4

p1 = posterior_difference_plot(
  recovery_noninvasive_normal, # Data object

  ## start using quotes
  "#A65041", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # Spacing between ticks in X axis
) 

p2 = posterior_cumulative_plot(
  recovery_noninvasive_normal, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  1, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial and previous evidence",
  subtitle = "Non-invasive ventilation subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

Invasive mechanical ventilation

priors_invasive_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "invasive ventilation"
  )

priors_invasive =
  rma(
    # Outcome
    measure = "RD", # risk difference
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = priors_invasive_raw %>% filter(outcome == "mortality"),

    slab = paste(trial),
    method = "REML"
  )

forest(priors_invasive)

set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size

data_prior_posterior_plot(
  recovery_invasive_normal, # Output from normal_approximation()
  
  # start using quotes
           "#364D55", # Color for the prior
           "#7EBAC2", # Color for RECOVERY
           "#A65041", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution - RECOVERY and previous evidence", # Title
           "Invasive mechanical ventilation subgroup", # Subtitle
           # stop using quotes
           -15, # X axis inferior limit
           10, # X axis superior limit
           5 # Spacing between ticks in X axis
  )
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlabel = "\nRisk difference (%)"
xlim_inf = -12
xlim_sup = 8

p1 = posterior_difference_plot(
  recovery_invasive_normal, # Data object

  ## start using quotes
  "#A65041", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # Spacing between ticks in X axis
) 

p2 = posterior_cumulative_plot(
  recovery_invasive_normal, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2, # Spacing between ticks in X axis
  "mortality" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial and previous evidence",
  subtitle = "Invasive mechanical ventilation subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is less than or equal to the effect size on the X‐axis.

Summary

### Create variables for the summary plot

set.seed = 123 # set seed for reproducibility (rnorm())
n = 10e4 # sampling size


## Non-informative priors

delta_beta_recovery_simple_oxygen =
  beta_recovery_simple_oxygen %>%
  # Calculate difference between groups
  summarise("Simple oxygen only" = toci - control)

delta_beta_recovery_non_invasive =
  beta_recovery_non_invasive %>% 
  summarise("Non-invasive ventilation" = toci - control)

delta_beta_recovery_invasive =
  beta_recovery_invasive %>% 
  summarise("Invasive mechanical ventilation" = toci - control)

# Bind all 3 tibbles together
delta_noninformative = 
  delta_beta_recovery_simple_oxygen %>% 
  bind_cols(delta_beta_recovery_non_invasive, delta_beta_recovery_invasive)

## Evidence-based priors

delta_recovery_simple_oxygen_normal =
  recovery_simple_oxygen_normal %>%
  summarise("Simple oxygen only" = rnorm(n,
                                         mean = post.mean,
                                         sd = post.sd)
            )

delta_recovery_noninvasive_normal =
  recovery_noninvasive_normal %>%
  summarise("Non-invasive ventilation" = rnorm(n,
                                               mean = post.mean,
                                               sd = post.sd)
            )

delta_recovery_invasive_normal =
  recovery_invasive_normal %>%
  summarise("Invasive mechanical ventilation" = rnorm(n,
                                              mean = post.mean,
                                              sd = post.sd)
            )

# Bind all 3 tibbles together
delta_evidence = 
  delta_recovery_simple_oxygen_normal %>% 
  bind_cols(delta_recovery_noninvasive_normal, delta_recovery_invasive_normal)
# Non-informative plot

# X axis
xlim_inf = -10
xlim_sup = 6
ticks_spacing = 2

# Assure subgroup order
# https://stackoverflow.com/questions/12774210/
# how-do-you-specifically-order-ggplot2-x-axis-instead-of-alphabetical-order

level_order = c("Invasive mechanical ventilation", "Non-invasive ventilation", "Simple oxygen only")

p1 = delta_noninformative %>%
  
  # make it tidy for ggplot
  pivot_longer(
    cols = c("Simple oxygen only":"Invasive mechanical ventilation"),
    names_to = "dist", values_to = "samples"
  ) %>%
  
  ggplot(aes(x = 100 * samples, y = factor(dist, level = level_order))) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x < 0)),
               .width = c(0.5, 0.95)
  ) +
  scale_fill_manual(values = c("gray80", "#7EBAC2")) +
  labs(
    x = "",
    y = "",
    title = "Non-informative priors\n"
  ) +
  scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    plot.title.position = "plot",
    legend.position = "none",
    plot.margin = margin(20, 25, 0, 20)
  ) +
  coord_cartesian(xlim = c(xlim_inf, xlim_sup))
# Evidence-based plot

p2 = delta_evidence %>%
  
  # make it tidy for ggplot
  pivot_longer(
    cols = c("Simple oxygen only":"Invasive mechanical ventilation"),
    names_to = "dist", values_to = "samples"
  ) %>%
  
  ggplot(aes(x = 100 * samples, y = factor(dist, level = level_order))) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x < 0)),
               .width = c(0.5, 0.95)
  ) +
  scale_fill_manual(values = c("gray80", "#A65041")) +
  labs(x = "\nRisk difference (%)",
       y = "",
       title = "Evidence-based priors\n"
  ) +
  scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(0, 25, 10, 20)
  ) +
  coord_cartesian(xlim = c(xlim_inf, xlim_sup))
p1 / p2 + plot_annotation(title = "Posterior distributions on mortality",
                          theme = theme(plot.title = element_text(size = 20)))
Interval bars depict 50% and 95% credible intervals.

Interval bars depict 50% and 95% credible intervals.

Hospital discharge

Regarding this outcome, negative risk differences mean benefit.


Using a non-informative prior

Simple oxygen only

### Generate posterior beta distributions using a noninformative beta(1,1) prior

set.seed = 123 # set seed for reproducibility (rbeta() function)

n = 10e4 # To determine sampling size

# Here I use the uninformative_posterior_beta() in which you can find in the
# following path: 
# here("final_analyses", "script", "functions", "beta_distributions.R")

beta_recovery_simple_oxygen = uniformative_posterior_beta( 
  recovery_simple_oxygen, # data
  "discharge" # outcome
) 

### Plot!

# I will use my own functions to standardize my plots

risk_comparison_plot(
  beta_recovery_simple_oxygen, # Data object

  ### start using quotes
  
  "gray50", # Color code for control group
  "#5C87C3", # Color code for other group
  "\nRisk (%)", # X axis label
  "RECOVERY trial", # Title
  "Simple oxygen only subgroup\n", # Subtitle
  
  ### stop using quotes

  58 , # X axis inferior limit
  76, # X axis superior limit
  2 # X axis spacing for ticks
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlabel = "\nRisk difference (%)"
xlim_inf = -2
xlim_sup = 14

p1 = risk_difference_plot(
  beta_recovery_simple_oxygen, # Data object

  ## start using quotes
  "#9D95C6", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # X axis label
) 

p2 = cumulative_risk_difference_plot(
  beta_recovery_simple_oxygen, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "RECOVERY trial",
  subtitle = "Risk difference between groups on simple oxygen only"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

Non-invasive ventilation

### Generate posterior beta distributions using a noninformative beta(1,1) prior

set.seed = 123 # set seed for reproducibility (rbeta() function)

n = 10e4 # To determine sampling size

# Here I use the uninformative_posterior_beta() in which you can find in the
# following path: 
# here("final_analyses", "script", "functions", "beta_distributions.R")

beta_recovery_non_invasive = uniformative_posterior_beta( 
  recovery_non_invasive, # data
  "discharge" # outcome
) 

### Plot!

# I will be using my own functions to standardize my plots

risk_comparison_plot(
  beta_recovery_non_invasive, # Data object

  ### start using quotes
  
  "gray50", # Color code for control group
  "#5C87C3", # Color code for other group
  "\nRisk (%)", # X axis label
  "RECOVERY trial", # Title
  "Non-invasive ventilation subgroup\n", # Subtitle
  
  ### stop using quotes

  34, # X axis inferior limit
  52, # X axis superior limit
  2 # Spacing between ticks in X axis
)
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlabel = "\nRisk difference (%)"
xlim_inf = -2
xlim_sup = 16

p1 = risk_difference_plot(
  beta_recovery_non_invasive, # Data object

  ## start using quotes
  "#9D95C6", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # Spacing between ticks in X axis
) 

p2 = cumulative_risk_difference_plot(
  beta_recovery_non_invasive, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "RECOVERY trial",
  subtitle = "Risk difference between groups on non-invasive ventilation"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

Invasive mechanical ventilation

### Generate posterior beta distributions using a noninformative beta(1,1) prior

set.seed = 123 # set seed for reproducibility (rbeta() function)

n = 10e4 # To determine sampling size

# Here I use the uninformative_posterior_beta() in which you can find in the
# following path: 
# here("final_analyses", "script", "functions", "beta_distributions.R")

beta_recovery_invasive = uniformative_posterior_beta( 
  recovery_invasive, # data
  "discharge" # outcome
) 

### Plot!

# I will be using my own functions to standardize my plots

risk_comparison_plot(
  beta_recovery_invasive, # Data object

  ### start using quotes
  
  "gray50", # Color code for control group
  "#5C87C3", # Color code for other group
  "\nRisk (%)", # X axis label
  "RECOVERY trial", # Title
  "Invasive mechanical ventilation subgroup\n", # Subtitle
  
  ### stop using quotes

  10, # X axis inferior limit
  26, # X axis superior limit
  2 # Spacing between ticks in X axis
) + 
  
  # Extra function to better limits the axis in this case
  coord_cartesian(c(10, 26))
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlabel = "\nRisk difference (%)"
xlim_inf = -8
xlim_sup = 12

p1 = risk_difference_plot(
  beta_recovery_invasive, # Data object

  ## start using quotes
  "#9D95C6", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  4 # Spacing between ticks in X axis
) 

p2 = cumulative_risk_difference_plot(
  beta_recovery_invasive, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "RECOVERY trial",
  subtitle = "Risk difference between groups on invasive mechanical ventilation"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

Using evidence-based priors

Simple oxygen only

priors_simple_oxygen_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "simple oxygen"
  )

priors_simple_oxygen =
  rma(
    # Outcome
    measure = "RD", # risk difference
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = priors_simple_oxygen_raw %>% filter(outcome == "discharge"),

    slab = paste(trial),
    method = "REML"
  )

forest(priors_simple_oxygen)

set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size

xlabel = "\nRisk difference (%)"

data_prior_posterior_plot(
  recovery_simple_oxygen_normal, # Output from normal_approximation()
  
  # start using quotes
           "#364D55", # Color for the prior
           "#9D95C6", # Color for RECOVERY
           "#DBA237", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution: RECOVERY and previous evidence", # Title
           "Simple oxygen only subgroup", # Subtitle
           # stop using quotes
           -10, # X axis inferior limit
           20, # X axis superior limit
           5 # Spacing between ticks in X axis
  )
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlim_inf = 0
xlim_sup = 14

p1 = posterior_difference_plot(
  recovery_simple_oxygen_normal, # Data object

  ## start using quotes
  "#DBA237", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # Spacing between ticks in X axis
) 

p2 = posterior_cumulative_plot(
  recovery_simple_oxygen_normal, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial and previous evidence",
  subtitle = "Simple oxygen only subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

Non-invasive ventilation

priors_noninvasive_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "non-invasive ventilation"
  )

priors_noninvasive =
  rma(
    # Outcome
    measure = "RD", # risk difference
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = priors_noninvasive_raw %>% filter(outcome == "discharge"),

    slab = paste(trial),
    method = "REML"
  )

forest(priors_noninvasive)

set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size

data_prior_posterior_plot(
  recovery_noninvasive_normal, # Output from normal_approximation()
  
  # start using quotes
           "#364D55", # Color for the prior
           "#9D95C6", # Color for RECOVERY
           "#DBA237", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution - RECOVERY and previous evidence", # Title
           "Non-invasive ventilation subgroup", # Subtitle
           # stop using quotes
           -10, # X axis inferior limit
           15, # X axis superior limit
           5 # Spacing between ticks in X axis
  )
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlim_inf = -2
xlim_sup = 14

p1 = posterior_difference_plot(
  recovery_noninvasive_normal, # Data object

  ## start using quotes
  "#DBA237", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2 # Spacing between ticks in X axis
) 

p2 = posterior_cumulative_plot(
  recovery_noninvasive_normal, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  1, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial and previous evidence",
  subtitle = "Non-invasive ventilation subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

Invasive mechanical ventilation

priors_invasive_raw =
  d %>%
  filter(
    trial != "RECOVERY",
    oxygen == "invasive ventilation"
  )

priors_invasive =
  rma(
    # Outcome
    measure = "RD", # risk difference
    
    # Tocilizumab group
    ai = events_toci,
    n1i = total_toci,

    # Control group
    ci = events_control,
    n2i = total_control,
    
    data = priors_invasive_raw %>% filter(outcome == "discharge"),

    slab = paste(trial),
    method = "REML"
  )

forest(priors_invasive)

set.seed = 123 # set seed for reproducibility (rnorm() function)
n = 10e4 # sampling size

data_prior_posterior_plot(
  recovery_invasive_normal, # Output from normal_approximation()
  
  # start using quotes
           "#364D55", # Color for the prior
           "#9D95C6", # Color for RECOVERY
           "#DBA237", # Color for the posterior
           xlabel, # X axis label
           "Posterior distribution - RECOVERY and previous evidence", # Title
           "Invasive mechanical ventilation subgroup", # Subtitle
           # stop using quotes
           -10, # X axis inferior limit
           25, # X axis superior limit
           5 # Spacing between ticks in X axis
  )
While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI.

While the darker shaded area depicts the 50% credible interval (CI), the interval bar depicts the 95% CI. 

xlabel = "\nRisk difference (%)"
xlim_inf = -8
xlim_sup = 12

p1 = posterior_difference_plot(
  recovery_invasive_normal, # Data object

  ## start using quotes
  "#DBA237", # fill color code
  xlabel, # X axis label
  ### stop using quotes

  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  4 # Spacing between ticks in X axis
) 

p2 = posterior_cumulative_plot(
  recovery_invasive_normal, # Data object
  xlabel, # X axis label (in quotes)
  xlim_inf, # X axis inferior limit
  xlim_sup, # X axis superior limit
  2, # Spacing between ticks in X axis
  "discharge" # Outcome (within quotes)
  )
p1 + p2 + plot_annotation(
  title = "Posterior distribution: RECOVERY trial and previous evidence",
  subtitle = "Invasive mechanical ventilation subgroup"
)
The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

The interval bar depicts the 95% credible interval. The cumulative posterior distribution corresponds to the probabilities that the risk difference (RD) is greater than or equal to the effect size on the X‐axis.

Summary

### Create variables for the summary plot

set.seed = 123 # set seed for reproducibility (rnorm())
n = 10e4 # sampling size


## Non-informative priors

delta_beta_recovery_simple_oxygen =
  beta_recovery_simple_oxygen %>%
  # Calculate difference between groups
  summarise("Simple oxygen only" = toci - control)

delta_beta_recovery_non_invasive =
  beta_recovery_non_invasive %>% 
  summarise("Non-invasive ventilation" = toci - control)

delta_beta_recovery_invasive =
  beta_recovery_invasive %>% 
  summarise("Invasive mechanical ventilation" = toci - control)

# Bind all 3 tibbles together
delta_noninformative = 
  delta_beta_recovery_simple_oxygen %>% 
  bind_cols(delta_beta_recovery_non_invasive, delta_beta_recovery_invasive)

## Evidence-based priors

delta_recovery_simple_oxygen_normal =
  recovery_simple_oxygen_normal %>%
  summarise("Simple oxygen only" = rnorm(n,
                                         mean = post.mean,
                                         sd = post.sd)
            )

delta_recovery_noninvasive_normal =
  recovery_noninvasive_normal %>%
  summarise("Non-invasive ventilation" = rnorm(n,
                                               mean = post.mean,
                                               sd = post.sd)
            )

delta_recovery_invasive_normal =
  recovery_invasive_normal %>%
  summarise("Invasive mechanical ventilation" = rnorm(n,
                                                      mean = post.mean,
                                                      sd = post.sd)
            )

# Bind all 3 tibbles together
delta_evidence = 
  delta_recovery_simple_oxygen_normal %>% 
  bind_cols(delta_recovery_noninvasive_normal, delta_recovery_invasive_normal)
# Non-informative plot

# X axis
xlim_inf = -8
xlim_sup = 14
ticks_spacing = 2

# Assure subgroup order
# https://stackoverflow.com/questions/12774210/
# how-do-you-specifically-order-ggplot2-x-axis-instead-of-alphabetical-order

level_order = c("Invasive mechanical ventilation", "Non-invasive ventilation", "Simple oxygen only")

p1 = delta_noninformative %>%
  
  # make it tidy for ggplot
  pivot_longer(
    cols = c("Simple oxygen only":"Invasive mechanical ventilation"),
    names_to = "dist", values_to = "samples"
  ) %>%
  
  ggplot(aes(x = 100 * samples, y = factor(dist, level = level_order))) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x > 0)),
               .width = c(0.5, 0.95)
  ) +
  scale_fill_manual(values = c("gray80", "#9D95C6")) +
  labs(
    x = "",
    y = "",
    title = "Non-informative priors\n"
  ) +
  scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    plot.title.position = "plot",
    legend.position = "none",
    plot.margin = margin(20, 25, 0, 20)
  ) +
  coord_cartesian(xlim = c(xlim_inf, xlim_sup))
# Evidence-based plot

p2 = delta_evidence %>%
  
  # make it tidy for ggplot
  pivot_longer(
    cols = c("Simple oxygen only":"Invasive mechanical ventilation"),
    names_to = "dist", values_to = "samples"
  ) %>%
  
  ggplot(aes(x = 100 * samples, y = factor(dist, level = level_order))) +

  # https://mjskay.github.io/ggdist/articles/slabinterval.html
  stat_halfeye(aes(fill = stat(x > 0)),
               .width = c(0.5, 0.95)
  ) +
  scale_fill_manual(values = c("gray80", "#DBA237")) +
  labs(x = "\nRisk difference (%)",
       y = "",
       title = "Evidence-based priors\n"
  ) +
  scale_x_continuous(breaks = seq(from = xlim_inf, to = xlim_sup, ticks_spacing)) +
  scale_y_discrete(expand = c(0, 0.3)) +
  theme(
    axis.ticks.x = element_blank(),
    axis.ticks.y = element_blank(),
    panel.background = element_blank(),
    panel.grid.major.x = element_line(color = "gray80", size = 0.3),
    legend.position = "none",
    plot.title.position = 'plot',
    plot.margin = margin(0, 25, 10, 20)
  ) +
  coord_cartesian(xlim = c(xlim_inf, xlim_sup))
p1 / p2 + plot_annotation(title = "Posterior distributions on hospital discharge",
                          theme = theme(plot.title = element_text(size = 20)))
Interval bars depict 50% and 95% credible intervals.

Interval bars depict 50% and 95% credible intervals. 



sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] metafor_2.4-0   Matrix_1.3-2    here_1.0.1      patchwork_1.1.1
##  [5] ggdist_2.4.0    tidybayes_2.3.1 flextable_0.6.4 forcats_0.5.1  
##  [9] stringr_1.4.0   dplyr_1.0.5     purrr_0.3.4     readr_1.4.0    
## [13] tidyr_1.1.3     tibble_3.1.1    ggplot2_3.3.3   tidyverse_1.3.1
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-152          fs_1.5.0              lubridate_1.7.10     
##  [4] httr_1.4.2            rprojroot_2.0.2       tools_4.0.4          
##  [7] backports_1.2.1       bslib_0.2.4           utf8_1.2.1           
## [10] R6_2.5.0              DBI_1.1.1             colorspace_2.0-0     
## [13] withr_2.4.1           tidyselect_1.1.0      curl_4.3             
## [16] compiler_4.0.4        cli_2.4.0             rvest_1.0.0          
## [19] arrayhelpers_1.1-0    xml2_1.3.2            officer_0.3.17       
## [22] sass_0.3.1            scales_1.1.1          systemfonts_1.0.1    
## [25] digest_0.6.27         rmarkdown_2.7         base64enc_0.1-3      
## [28] pkgconfig_2.0.3       htmltools_0.5.1.1     highr_0.9            
## [31] dbplyr_2.1.1          fastmap_1.1.0         rlang_0.4.10         
## [34] readxl_1.3.1          rstudioapi_0.13       jquerylib_0.1.3      
## [37] generics_0.1.0        farver_2.1.0          svUnit_1.0.3         
## [40] jsonlite_1.7.2        zip_2.1.1             distributional_0.2.2 
## [43] magrittr_2.0.1        Rcpp_1.0.6            munsell_0.5.0        
## [46] fansi_0.4.2           gdtools_0.2.3         lifecycle_1.0.0      
## [49] stringi_1.5.3         yaml_2.2.1            plyr_1.8.6           
## [52] grid_4.0.4            crayon_1.4.1          lattice_0.20-41      
## [55] haven_2.4.0           hms_1.0.0             knitr_1.32           
## [58] pillar_1.6.0          uuid_0.1-4            reprex_2.0.0         
## [61] glue_1.4.2            evaluate_0.14         V8_3.4.0             
## [64] data.table_1.14.0     modelr_0.1.8          vctrs_0.3.7          
## [67] cellranger_1.1.0      gtable_0.3.0          assertthat_0.2.1     
## [70] cachem_1.0.4          xfun_0.22             broom_0.7.6          
## [73] coda_0.19-4           memoise_2.0.0         ellipsis_0.3.1       
## [76] rdocsyntax_0.4.1.9000